1. show where your plots are

2. show how variables are distributed spatially

3. spatial analyses

4. Many, many tools

1. Workflow

2. Quite efficient

3. Interface

R for everything (almost)

1. import your data
2. format your data
3. analyze your data
4. visualize your data
5. export your data
6. create your own functions/packages

1. well-defined classes
   
2. many formats (+ convert different file)
   
3. operations on geometries
   
4. statistical analyses
   
5. visualization

1. C/C++: Geos, proj4, Gdal
2. JavaScript: Leaflet
3. otherGIS: RQGIS, rgrass7, RgoogleMaps
4. retrieve data: osmdata

1- Define classes

• sp
• raster
• sf

2- Topology operations

• rgeos
• sf

3- Import / export

• rgdal
• sf

4- Analyses

• spdep
• dismo

5- Creating and editing maps

• graphics
• ggplot
• interactive plot / web-based plot
• a blog post about this

6- But

• Geo-referencing
• Visualizing large spatial objects

library(sp)

## Loading required package: methods

mysp <- SpatialPointsDataFrame(
  coords = mydata[,1:2],
  data = mydata[,3:4],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

class(mysp)

## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"

isS4(mysp)

## [1] TRUE

slotNames(mysp)

## [1] "data"        "coords.nrs"  "coords"      "bbox"        "proj4string"

mysp@proj4string

## CRS arguments:
##  +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0

head(mysp@data)

##         var1     var2
## 1 -0.9424900 8.822271
## 2  0.4211175 9.058452
## 3  2.1195748 4.783351
## 4  0.3377636 8.061516
## 5  0.8181175 4.986830
## 6 -0.5898903 4.639604

(mysp2 <- spTransform(mysp, CRS=CRS("+proj=merc +ellps=GRS80")))

##            coordinates        var1      var2
## 1  (-8914219, 5318662) -0.94249004 8.8222714
## 2  (-9081556, 5327927)  0.42111746 9.0584523
## 3  (-9047620, 5200757)  2.11957478 4.7833508
## 4  (-9098812, 5330904)  0.33776361 8.0615156
## 5  (-9076630, 5326812)  0.81811746 4.9868304
## 6  (-8931971, 5229491) -0.58989032 4.6396036
## 7  (-9049053, 5341975) -0.14429222 7.9327105
## 8  (-9082845, 5365187)  0.73672200 4.9249723
## 9  (-9010743, 5223804)  0.24897363 4.0290142
## 10 (-9042984, 5408682) -0.08106143 3.0501115
## 11 (-8992699, 5358424) -0.46859518 0.6339195
## 12 (-8945265, 5158475)  0.70483256 1.5400158
## 13 (-8976821, 5431270) -0.62466071 0.4572709
## 14 (-9095845, 5201313)  1.93421826 0.9320450
## 15 (-9110967, 5255899)  0.77250278 6.2389721
## 16 (-9024592, 5159324) -1.09241480 7.7116822
## 17 (-8935065, 5326269) -0.41269577 4.4730041
## 18 (-9070013, 5264491) -0.32887855 4.7979019
## 19 (-9011305, 5231279) -1.14365387 1.1729072
## 20 (-9039393, 5432599)  1.54523563 7.8279072

0- SpatialGrid / SpatialPixel (package sp)

3- RasterLayer (package raster)

4- RasterStack (package raster)

5- RasterBrick (package raster)

library(raster)

## 
## Attaching package: 'raster'

## The following object is masked from 'package:magrittr':
## 
##     extract

ras1 <- raster(matrix(runif(100*100,0,10),ncol=100,nrow=100),
    crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"),
    xmn=-82, xmx=+80, ymn=+42, ymx=+44)

ras1

## class       : RasterLayer 
## dimensions  : 100, 100, 10000  (nrow, ncol, ncell)
## resolution  : 1.62, 0.02  (x, y)
## extent      : -82, 80, 42, 44  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : layer 
## values      : 0.0001147902, 9.997746  (min, max)

head(getData("ISO3"))

##   ISO3        NAME
## 1  ABW       Aruba
## 2  AFG Afghanistan
## 3  AGO      Angola
## 4  AIA    Anguilla
## 5  ALA       Åland
## 6  ALB     Albania

## Country level:
mapBEL0 <- getData(name="GADM", country="BEL", path="./assets", level=0)
mapBEL1 <- getData(name="GADM", country="BEL", path="./assets", level=1)
tminW <- getData(name="worldclim", var="tmin", res=10, path="./assets")
mapBEL0

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 2.555356, 6.40787, 49.49722, 51.50382  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 68
## names       : OBJECTID, ID_0, ISO, NAME_ENGLISH, NAME_ISO, NAME_FAO,              NAME_LOCAL, NAME_OBSOLETE,            NAME_VARIANTS, NAME_NONLATIN, NAME_FRENCH, NAME_SPANISH, NAME_RUSSIAN, NAME_ARABIC, NAME_CHINESE, ... 
## value       :        1,   23, BEL,      Belgium,  BELGIUM,  Belgium, Belgique|Belgie|Belgien,              , Belgique|Belgien|Belgium,              ,   Belgique ,     Bélgica ,     Бельгия ,      بلجيكا,      比利时 , ...

library(rgdal)

## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: /usr/share/gdal/2.1
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.2-5

1. writeOGR / writeGDAL: to write spatial objects
2. readOGR/readGDAL: to read spatial files

head(ogrDrivers())

##         name                     long_name write  copy isVector
## 1 AeronavFAA                   Aeronav FAA FALSE FALSE     TRUE
## 2 AmigoCloud                    AmigoCloud  TRUE FALSE     TRUE
## 3     ARCGEN             Arc/Info Generate FALSE FALSE     TRUE
## 4     AVCBin      Arc/Info Binary Coverage FALSE FALSE     TRUE
## 5     AVCE00 Arc/Info E00 (ASCII) Coverage FALSE FALSE     TRUE
## 6        BNA                     Atlas BNA  TRUE FALSE     TRUE

writeOGR(mysp, dsn="./assets", layer="mypoints",
    driver="ESRI Shapefile", overwrite_layer=TRUE)

mysp2 <- readOGR(dsn="assets/", layer="mypoints")

## OGR data source with driver: ESRI Shapefile 
## Source: "assets/", layer: "mypoints"
## with 20 features
## It has 2 fields

## Roads find on http://www.diva-gis.org/Data
canroads <- readOGR(dsn="assets/", layer="CAN_roads")

## OGR data source with driver: ESRI Shapefile 
## Source: "assets/", layer: "CAN_roads"
## with 16915 features
## It has 5 fields

Note: there are also functions to import/export raster in raster

library(rgeos)

## rgeos version: 0.3-26, (SVN revision 560)
##  GEOS runtime version: 3.5.1-CAPI-1.9.1 r4246 
##  Linking to sp version: 1.2-5 
##  Polygon checking: TRUE

buf <- gBuffer(mapBEL0, width=0.5)

## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates

diff <- gDifference(gBuffer(mapBEL0, width=0.5), gBuffer(mapBEL0, width=0.1))

## Warning in gBuffer(mapBEL0, width = 0.5): Spatial object is not projected;
## GEOS expects planar coordinates

## Warning in gBuffer(mapBEL0, width = 0.1): Spatial object is not projected;
## GEOS expects planar coordinates

library(mapview)

NB: it uses the leaflet package.

mapview(mysp, cex = 'var2')@map

mapview(mapBEL1)@map

• https://insileco.github.io/Elections_fr-042017/

• https://insileco.github.io/fr-departement-election/

plot(mapBEL0)

class(tminW)

## [1] "RasterStack"
## attr(,"package")
## [1] "raster"

plot(tminW)

plot(mapBEL0, border='grey15', col='#E6E6E6', lwd=1.6)
plot(mapBEL1, lty=2, lwd=0.9, add=T)
points(4.3513, 50.8471, pch=19, col="#27df9d")
text(4.3513, 50.8471, text="Brussel", pos=3)

Simple features or simple feature access refers to a formal standard (ISO 19125-1:2004) that describes how objects in the real world can be represented in computers, with emphasis on the spatial geometry of these objects

https://cran.r-project.org/web/packages/sf/vignettes/sf1.html

1- Have a look at https://github.com/r-spatial/sf

2- Install R packages

install.packages(c("sf", "sp", "rgeos", "rgdal", "raster", "mapview"))

• the .csv file below is available in the repository

lakedf <- read.csv('assets/lakeOnt.csv')
head(lakedf)

##   X LAKEID LONGITUDE_DD LATITUDE_DD Area_ha
## 1 1 L32173    -80.03361    49.90417   174.9
## 2 2 L21001    -91.68278    48.68028   130.2
## 3 3 L44001    -84.19667    46.91861   279.6
## 4 4 L25002    -90.43111    48.16639   114.1
## 5 5 L44002    -84.30528    47.06000   137.3
## 6 6 L21003    -91.34528    48.24806  2993.9

library(sp)
lakesp <- SpatialPointsDataFrame(
  coords = lakedf[,3:4],
  data = lakedf[,c(1,4)],
  proj4string = CRS(
    "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"
  )
)

lakesp

## class       : SpatialPointsDataFrame 
## features    : 546 
## extent      : -95.06111, -76.01194, 43.05722, 54.52056  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0 
## variables   : 2
## names       :   X, LATITUDE_DD 
## min values  :   1,    43.05722 
## max values  : 546,    54.52056

altCAN <- getData(name="alt", country="CAN", path="./assets/")
altCAN

## class       : RasterLayer 
## dimensions  : 4992, 10620, 53015040  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : -141.1, -52.6, 41.6, 83.2  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
## data source : /home/kevcaz/Github/Tutorials/mapsWithR/docs/assets/CAN_msk_alt.grd 
## names       : CAN_msk_alt 
## values      : -108, 5800  (min, max)

bouCAN <- getData(country='CAN', level=1, path="./assets/")
bouCAN

## class       : SpatialPolygonsDataFrame 
## features    : 15 
## extent      : -141.0069, -52.61889, 41.67693, 83.11042  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : OBJECTID, ID_0, ISO, NAME_0, ID_1,  NAME_1, HASC_1, CCN_1, CCA_1,     TYPE_1, ENGTYPE_1, NL_NAME_1,                                       VARNAME_1 
## min values  :        1,   42, CAN, Canada,    1, Alberta,  CA.AB,    NA,    10,   Province,  Province,          ,                                                 
## max values  :       15,   42, CAN, Canada,   13,   Yukon,  CA.YT,    NA,    62, Territoire, Territory,          , Yukon Territory|Territoire du Yukon|Yukon|Yuk¢n

bouCAN@data[,1:8]

##    OBJECTID ID_0 ISO NAME_0 ID_1                    NAME_1 HASC_1 CCN_1
## 1         1   42 CAN Canada    1                   Alberta  CA.AB    NA
## 2         2   42 CAN Canada    2          British Columbia  CA.BC    NA
## 3         3   42 CAN Canada    3                  Manitoba  CA.MB    NA
## 4         4   42 CAN Canada    4             New Brunswick  CA.NB    NA
## 5         5   42 CAN Canada    5 Newfoundland and Labrador  CA.NF    NA
## 6         6   42 CAN Canada    6     Northwest Territories  CA.NT    NA
## 7         7   42 CAN Canada    7               Nova Scotia  CA.NS    NA
## 8         8   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 9         9   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 10       10   42 CAN Canada    8                   Nunavut  CA.NU    NA
## 11       11   42 CAN Canada    9                   Ontario  CA.ON    NA
## 12       12   42 CAN Canada   10      Prince Edward Island  CA.PE    NA
## 13       13   42 CAN Canada   11                    Québec  CA.QC    NA
## 14       14   42 CAN Canada   12              Saskatchewan  CA.SK    NA
## 15       15   42 CAN Canada   13                     Yukon  CA.YT    NA

bouCAN[11,]

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -95.15598, -74.34605, 41.67693, 56.86972  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 13
## names       : OBJECTID, ID_0, ISO, NAME_0, ID_1,  NAME_1, HASC_1, CCN_1, CCA_1,   TYPE_1, ENGTYPE_1, NL_NAME_1,    VARNAME_1 
## value       :       11,   42, CAN, Canada,    9, Ontario,  CA.ON,    NA,    35, Province,  Province,          , Upper Canada

plot(crop(altCAN, bouCAN[11,]@bbox))
par(mar=c(0,0,0,0))

altONT <- rasterize(bouCAN[11,], crop(altCAN, bouCAN[11,]@bbox), mask=TRUE)
mapview(altONT)+mapview(lakesp)

par(mar=c(0,0,0,0))
plot(altONT)
plot(lakesp, add=T)

• task view: https://cran.r-project.org/web/views/Spatial.html; also following to this https://cran.r-project.org/web/views/

install.packages("ctv")
install.views("Spatial")

• Document par Nicolas Casajus
• http://r-spatial.org/

• free data at country level
• Quebec free data
• find more spatial data
• create shapefile on line